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ABSTRACT 

We use the recently completed one billion particle Via Lactea II ACDM simulation to 
investigate local properties like density, mean velocity, velocity dispersion, anisotropy, 
orientation and shape of the velocity dispersion ellipsoid, as well as structure in ve- 
locity space of dark matter haloes. We show that at the same radial distance from the 
halo centre, these properties can deviate by orders of magnitude from the canonical, 
spherically averaged values, a variation that can only be partly explained by triaxiality 
and the presence of subhaloes. The mass density appears smooth in the central relaxed 
regions but spans four orders of magnitude in the outskirts, both because of the pres- 
ence of subhaloes as well as of underdense regions and holes in the matter distribution. 
In the inner regions the local velocity dispersion ellipsoid is aligned with the shape 
ellipsoid of the halo. This is not true in the outer parts where the orientation becomes 
more isotropic. The clumpy structure in local velocity space of the outer halo can not 
be well described by a smooth multivariate normal distribution. Via Lactea II also 
shows the presence of cold streams made visible by their high 6D phase space density. 
Generally, the structure of dark matter haloes shows a high degree of graininess in 
phase space that cannot be described by a smooth distribution function. 

Key words: Galaxies: haloes — Galaxies: structure — Galaxies: kinematics and 
dynamics — dark matter — methods: TV-body simulations — methods: numerical 



1 INTRODUCTION 

Spherical averaging is a commonly used method to describe 
the characteristics of dark matter haloes that form by grav- 
itational collapse in a cosmological environment. For exam- 
ple one of the basic characteristics of a dark matter halo 
is its spherically-averaged density profile which can be well 
described by a smooth function within resolved scales (e.g. 



Navarro et al.|1996| [Moore et al.|1998| |Diemand et~aL[ 2005 , 



2007||2008| jStadel et al.||2008| ). Other examples include the 
velocity dispersion and anisotropy profiles. 

Consider a set of observers at a given distance from the 
halo centre and capable of measuring local properties of the 
halo (e.g., density, velocity dispersion, velocity anisotropy). 
In a smooth, spherically symmetric halo, these measure- 
ments would yield identical values, up to statistical fluc- 
tuations. However, it is well known that dark matter haloes 
that form in pure dark matter simulations, which neglect 
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possible baryonic effects in the centre, are in general triax- 
ial: close to prolate in the central part and becoming rounder 
in the outskirts of the halo (e.g. Katz|199lj|Dubinski fc Carl- 



berg 199l1 |Allgood et aI.|2006||Bett et al.|2007||Kuhlen et al. 



2007). Such a shape variation obviously leads to a significant 



variance in local properties compared to the spherically av- 
eraged value at a given radius flKnebe fe Wiefiner 12006) . In 
addition, haloes have a high level of subhaloes which also ef- 
fect results of local measurements. But local deviations from 
a smooth model go beyond these shape and subhalo driven 
variations, as we demonstrate in this paper. For example, 
overdensities are expected due to the presence of subhaloes, 
but we also find underdense regions, which are unexpected 
in a smooth triaxial background halo with subhaloes (see 
section 



3.2) 



One aim of this paper is to quantify the degree of vari- 
ation between local and spherically averaged values for key 
quantities such as the density and velocity dispersion. A 
better understanding of the local variations of these prop- 
erties is essential in order to obtain a better description of 
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the phase space structure of dark matter haloes and their 
formation process. 

Such local variations can have profound consequences. 
For example, in dynamical models of dark matter haloes and 
galaxies it is often assumed that most properties are just a 
function of radius (e.g. jDehnen &; McLa ughlin 2005 Baes 



fe van Hesel|2007| |Van Hese et al.||2008| ). From the work 



presented here it becomes clear that this is not a valid as- 
sumption for the anisotropy (see section 3.5 ) and a lot of in- 



formation is smeared out by spherically averaging. Models 
based on simple analytical forms of the distribution func- 
tion, which are also used for generating iV-body realisations 
of dark matter haloes a nd galaxies (e.g. |Kazantzidis et aL] 
2004| |Zemp et aL||2008| |Zhang fe Magorrian||2008| ), often 
have the property that the local bulk motion is zero and the 
velocity distribution function is similar to a smooth mul- 
tivariate normal distribution (if it is not even deliberately 
set to a multivariate normal distribution as in Hernquist 
(1993)). Both of these characteristics are not valid in the 



outskirts of dark matter haloes that form in a hierarchical 



cosmological context (see sections 3.3 and 3.7). As we see, 
most of these effects are not taken into account in current 
dynamical models of dark matter haloes and galaxies and 
they will be a challenge to be modelled with more realistic 
distribution functions. 

In addition, it is also important to quantify these varia- 
tions locally at the Earth's position within the Galaxy in 
order to interpret results from present and future direct 



dark matter detection experiments on Earth (e.g. |Moore 
eTH1|200T] IStiff et al.|[^00T] [Hehrn et al.| [2002] [Stiff gg 



Widrow 2003; Kamionkowski & Koushiappas 2008; Vogels- 



berger et al.||2008| |Fairbairn fe Schwetz||2 008). While small 
subhalos contribute significantly to the indirect detection 



signal (e.g. Diemand et al. 2006 Kuhlen et al. 2008) and 



infall caustics have no significant effect (IDiemand & Kuhlen 



2008 ) , it is worthwhile also to investigate whether additional 



lumpiness from tidal debris might enhance the annihilation 



rate (see section 3.2I 1 



Also recent results from surveys like SDSS and RAVE 
need a better understanding of the structure of the outer 
halo. There is a lot of known structure in the Galactic stellar 
halo like for example the Sagittarius stream, the Monoceros 



stream or the Virgo overdensity (e.g. Martmez-Delgado et al. 
2007||Cole et al.||2008| |Casetti-Dinescu et al.||2008| ). Obvi- 



ously, these features are baryonic, but stars are collisionless 
too, and it's probable that these features originate from in- 
falling dark matter dominated objects. 

In general, a more detailed picture of the phase space 
structure of dark matter haloes is needed in order to under- 
stand and model their properties that are set by the hierar- 
chical formation process. This work is a further important 
step towards that aim. 

Here we present a study of the distribution of these 
properties as a function of galactocentric distance of a Milky 
Way size dark matter halo with data from the Via Lactea II 
project ( Diemand et al.||2008 ). In section [2] we give a sum- 



1 After submitting this work an independent, analytic study of 
this question (Afsh ordi et al.|2 008) has appeared. Our simulation 
results agree that such an enhancement is rather small (see section 
!T2!. 



mary of the simulation data used in this study and present 
the results in section [3] We then discuss our results and 
present the conclusions in section [4] In appendix [A] we give 
a comparison of the results with lower resolution simulations 
and different definitions of locality. 



2 SIMULATION AND DATA 

For the analysis presented here, we used data from the Via 
Lactea II simulation ( |Diemand et al.|2008| ) which simulates 
the assembly of a Milky Way size cold dark matter halo 
within a ACDM universe. The initial conditions at a starting 
redshift of z — 104.3 consist of a 40 comoving Mpc periodic 
box and were generated with a modified, parallel version 



of GRAFIC2 flBertschinger||2001| ). We use the traditional 
method of refining a region of interest with a large number 
of particles and leaving the rest on lower resolution so that 
we correctly account for the large scale tidal forces ( |Katz| 
\h White||1993||Bertschinger|[200T| . In total, the simulation 
consists of more than 10 high resolution particles with a 
mass of m p = 4098 Mq and a gravitational softening length 
of 40 pc. We used the WMAP 3-year cosmological parame- 
ters |Spergel et aLl|2007|) with Q M ,o 0.238, Q A ,o 0.762 



The time evolution until redshift z = was per- 



formed with the parallel tree-code PKDGRAV2 (Stadel 



2001). PKDGRAV2 uses a fast multipole expansion tech- 



nique in order to calculate the forces with hexadecapole pre- 
cision and a time-stepping scheme that is based on the true 
dynamical time of the particles with an accuracy parameter 
of 77d = 0.06 ( [Zemp et al.|[2007| [Stadel fe Zemp|2008) . The 
simulation used close to 10 CPUh on the Jaguar Cray XT3 
super-computer at the Oak Ridge National Laborator^ 

Throughout the paper, we only use data from the z — 
snapshot. We present all quantities in physical units and in 
a coordinate system centred on the particle with the deep- 
est potential within the dark matter halo. A detailed dis- 
cussion of the Via Lactea II simulation at z = can be 
found in IDiemand et al. ( 2008 ) and we report here only 



the main characteristics. The radius r2oob, where the en- 
closed density is 200 times the mean matter density (back- 
ground), is r2oob = 402.1 kpc and contains a total mass 
of M 200 b = M(r 20 ob) = 1.917 x 10 12 M . The maximum 
circular velocity of the halo is v C max = 201.3 km s _1 and 
is reached at the radius r Vcmax = 59.83 kpc. We determine 
the shape of the halo at 20 kpc by diagonalising the shape 
tensojj S with elements 



(i) 



2 http://www.ornl.gov 

3 Often the tensor S is incorrectly denoted as the inertia tensor. 
However, the tensor S does not give the correct relation between 
the angular momentum vector L and the angular velocity vector u? 
which is given by L = I cv where Uj = Xlfc mfc [ r fc^' _ ( x k)i( x k) j] 
is the correct definition of the inertia tensor. Of course, the tensors 
I and S have the same eigenvectors since they differ only by a 
diagonal tensor with one eigenvalue of algebraic multiciply of 3. 
Though, the eigenvalues of S and I are obviously different and 
have entirely different meanings. 
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Table 1. Summary of shell properties. 
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through an iterative procedure that adapts to the local shape 
as described in |Katz| ( |1991| ) so that the summation over the 
particles only contains those within the triaxial ellipsoid. We 
then rotate the whole halo so that the long axis is the x-axis, 
the intermediate axis is the y-axis, and the short axis is the 
z-axis. The resulting axis ratios at 20 kpc are b/a — 0.6154 
and c/a — 0.5227 where a ^ b ^ c are the square roots of 
the eigenvalues of the shape tensor. Of course, a dark matter 
halo is not a rigid body and the axis ratios as well as their 
orientations can change with radius, with haloes typically 
becoming rounder in their outer part. The orientation of 
the principle axis is constant with radius (to within a few 
degrees) between 20 kpc and 400 kpc. 



3 LOCAL PROPERTIES 



3.1 Procedure 



In order to measure local properties of the dark matter halo, 
we chose the following procedure. We consider six distinct 
distances: 8, 25, 50, 100, 200 and 400 kpc from the halo 
centre. At each distance we randomly sample A sam pie = 10 4 
spheres with a radius r sp h(r) depending only on the distance 
and chosen to include of O(10 3 ) particles. Specifically we 
set r sp h(8 kpc) = 0.5 kpc and calculate the radii at other 
distances by 



' sph 



(r) 



A sph I 



kpc) 



p(8 kpc) 
p(r) 



(2) 



where p(r) is the measured spherically averaged density at 
radius r. The number of sampling spheres is large enough so 
that even at 400 kpc the whole sky, as seen from the centre 
of the galaxy, is covered more than once. A summary of the 
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characteristics of the shells at different distances from the 
galactic centre is given in Table [I] In order to save space, 
in most cases we only show the plots for 8, 25, 100 and 400 
kpc in the following analysis. 

In order to assess the influence of the triaxial shape of 
the halo, we also consider 3 sub-sets of spheres whose centres 
lie within a cone of 20° along the long (x), the intermediate 
(y) respectively the short (z) axis, i.e. spheres that fulfil the 
condition 



|e sp h ■ ek\ > cos(20°) with k — x^yoiz 



(3) 



where e sp h is the unit vector in galactocentric coordinates 
towards the centre of the sphere and ek is the unit vector 
along the x-, y- respectively z-axis. It is expected that each 
sub-set contains approximately 6% of the total number of 
spheres and the actual numbers lie within the statistical 
range. The slightly twisted shape of the dark matter halo has 
only a small influence on the sub-sets along the different axes 
at different distances r since the deviations of the eigenvector 
directions are only a few degrees which is small compared 
to the cone opening angle of 20°. 

Furthermore we consider the sub-sample of spheres that 
are affected by subhaloes. We take a group catalogue of the 
Via Lactea II simulation that contains all subhaloes with 
a peak circular velocity greater than 2 km s _1 , of which 
there are 26182 within 500 kpc. These groups are identified 
by a friends-of-friends method in phase space (6DFOF) as 



described in |Diemand et al. ( 2006| . By a density profile fit- 
ting algorithm, we determine a size r s h for each subhalo. A 
sphere with centre at location x sp h is then deemed affected 
by a subhalo if the criterion 



|»^sph »Csh,fc| ^ ^sph ~h ^sh,fc 



(4) 



is fulfilled for any subhalo k with radius r s h,k and centred at 
x s h,k, i-e. the subhalo k at least partially overlaps the sphere. 
We give the fraction of spheres that are affected by subhaloes 
m Table [Q This fraction increases with distance from the 
centre of the galaxy, reaching a value of approximately 28% 
at 400 kpc. The fraction of subhalo- affected spheres is a 
lower limit since more subhaloes would survive in higher 
resolution simulations at all distances from the centre of the 
galaxy. 

Throughout this paper we use the following notation 
scheme: quantities calculated by averaging over particles 
within a sphere are denoted by q, quantities calculated by 
averaging over particles within a shell from r — r sp h(r) to 
r + r sp h(r) are denoted by q, and averages over the ensemble 
of iVsampie spheres at radius r are denoted by (g), where we 
also might use indices in order to denote ensemble averages 
over only sub-sets of the spheres, for any given quantity q. 

3.2 Local densities 

In Fig.[l] we plot the probability density function of the local 
mean density p within the spheres in logarithmic scale. The 
local density is simply defined by the mass within a sphere 
divided by the volume of the sphere. In the central region, 
the structure of the halo looks pretty smooth in density. At 
8 kpc, we find a maximum total spread of a factor zb 2 in 
local density with a standard deviation of a(p)/p w 0.26. 
This is due to the triaxiality of halo since the underdense 
spheres lie along the z-axis and the overdense spheres along 



the x-axis. In the outer part of the halo, however, the smooth 
picture does not apply. For example at 200 kpc and 400 kpc, 
underdense regions are found away from the short axis and 
the distributions of the sub-samples along the other axes are 
broader than in the centre of the halo. The possible range in 
densities at 400 kpc is close to 4 orders of magnitude with 
a peak probability at around 0.5 p(400 kpc), i.e. the typi- 
cal value for the local densities is below the spherical aver- 
age value. Beyond 100 kpc, the standard deviation becomes 
larger than the mean value, e.g. at 400 kpc a(p)/p ~ 5.3. 

Naively, one would expect a decrease in the density vari- 
ation due to shape in the outskirts of the halo, since its shape 
becomes rounder further out. However, this is exceeded by 
the steeper fall off of the density profile in the outskirts of 
the dark matter halo. The actual density contrast between 
the mean local density of the x-axis sample with respect 
to the mean local density of the z-axis, C p = (p) x / (p) z , 
increases from C P (S kpc) = 2.353 to C p (400 kpc) = 14.72. 

It is also striking that there are no high density spikes 
from subhaloes in the central region. This is just an artefact 
of our rather conservative definition of locality (i.e. the large 
size of the spheres). There are only very few surviving and 
well resolved subhaloes within 8 kpc in Via Lactea II. The 
subhaloes that survive close to the centre are small and com- 
pact due to tidal mass loss so that with our still relatively 
large definition of locality, their contribution to the mass of 
the spheres is small. By decreasing the size of the spheres 
the extremes of the distribution become more populated (see 
also appendix [A|. With higher resolution, one expects more 
subhaloes to survive which are even more compact. But one 
would also decrease the size of the spheres that define local- 
ity in our scheme. Hence, it is not clear what the outcome 
of these two competing effects would be. From simple an- 
alytical models, one expects some contribution from these 
small haloes ( |Kamionkowski fc Kou shiappas 2008) though 
the details depend on their volume filling factor. 

While high peaks from subhaloes are expected, under- 
dense regions away from the short axis of the shape ellip- 
soid are surprising. In order to further investigate these un- 
derdensities, we repeated the measurement of local prop- 
erties with spheres of 4 times smaller radii, i.e. r sp h/4. In 
these spheres, one would expect on average approximately 
21 particles in the inner region. We then simply count the 
small spheres that contain no particles and calculate from 
that the fraction of empty spheres (given in Table [TJ . As- 
suming balls-in-bins statistics (see appendix [b] for more de- 
tails), the probability of getting empty spheres is given by 



Pempty(2l) e 



7.583 x 10" 



Hence one would not 



expect to find any sphere of size r sp h/4 that is completely 
empty. The fact that approximately 2% of all these smaller 
spheres at 400 kpc are empty shows clearly that the outskirts 
of dark matter haloes are far from smooth in position space. 
Of course, to some degree this is due to the triaxial shape as 
actually most of the empty spheres in the outskirts are in low 
density regions along the z-axis. Taking the lower density 
into account, one would expect approximately 6 particles 
per sphere in the z-axis sample at 400 kpc. The probability 
of getting empty spheres IS Pempty (6)= e - 6 = 2.479 xl0- 3 . 
The measured fraction of empty spheres in the z-sample is 
/empty, z = 6.656 x 10 -2 which is still clearly much higher. 

We quantify this further with a statistical indicator that 
measures inequality: the Gini coefficient. The Gini coeffi- 
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Figure 1. Probability density functions of the local density p at different galactocentric distances r normalised by the spherically averaged 
value p (solid line/grey). Additionally, we plot the sub-samples along the different axes of the triaxial halo: x-axis sample (dash-dotted 
line/blue), y-axis sample (long-dashed line/red) and z-axis sample (short-dashed line/orange). We also plot the subhalo-affected sample 
(long-short-dashed line/green). The ensemble average (p) over the spheres is marked with a solid black square and the standard deviation 
range (p) ± cr(p) is marked with a solid black hexagon. The shell value p is marked with an open black circle. 



cient is defined as G = 1 — 2 A where A is the area under the 



1905 



Lorenz curve of the probability distribution (Lorenz 
Gini|1912 ). For the discrete probability density functions of 
the local mean density this results in 



G(p) 



with 



(5) 



^sample 



60 = and 

Pi(p) = Pi(p)Alog(p/p) . 



(6) 



(7) 



with Alog(p/p) being the bin width in logarithmic scale. 
The range of the Gini coefficient is between and 1. A Gini 
coefficient of means that all the spheres have equal density 
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whereas a Gini coefficient of 1 means that all the mass is in 
one sphere. As can be seen in Table ^ the Gini coefficient 
increases with galactocentric distance showing again that 
the outskirts of dark matter haloes have a clumpy structure. 

This clumpy structure is due to the hierarchical build 
up of dark matter haloes by accretion of subhaloes. These 
subhaloes leave tidal streams along their orbits when they 
are falling into the host halo. In the outer parts of the halo 
we find many overlapping streams locally (see also section 
3.7) that originate from the many subhaloes that passed 
a certain region. Combined with the long local dynamical 
time-scale in the outskirts of haloes which inhibits effective 
mixing, this leads to the clumpy structure and large spread 
in density we see on small local scales in the outer parts of 
dark matter haloes. 

The tidal debris can have an effect on the local dark 
matter annihilation (Afsh ordi et al.|2008| ). We can quantify 
this by the boost factor B = (p z ) / {pY ■ We can only say 
something about local boost factors since by averaging over 
a spherical shell, the boost factor would be mainly driven by 
the triaxiality. The biggest effect is expected in the outskirts 
of the halo since there we find the highest degree of dumpi- 
ness. In order just to measure the influence of the debris 
structure and exclude the effect of the subhaloes, we only 
include spheres that are not subhalo- affected and get values 
from ^,^(400 kpc) = 1. 319 to £ x , noa/t (4 kpc ) = 3.196. 
This is in agreement with Afshordi et al. (2008) who get 



of 0(1) for the boost factor originating in the structure 
of the debris. If we include the subhalo- affected sample 
we get values in the range from B y (400 kpc) = 6.046 to 
B x (400 kpc) = 16.15 locally, showing showing that even in 
the outer halo most of the dumpiness comes from the subha- 
los themselves and not from their more diffuse tidal debris. 

Generally, all the probability density function plots in 
this paper are lower limits for a possible real variation of 
local properties in nature since numerical effects due to nu- 
merical under-resolving (even in Via Lactea II) lead to ar- 
tificial heating, smoothing, and in general loss of structure. 
A more detailed discussion of numerical effects follows in 
appendix [A] 

3.3 Local mean velocities 

In this section we investigate the properties of the local mean 
velocities, which we calculate by number weighted averages 
within the spheres. For this purpose, we split the velocity 
vector field in a radial, azimuthal (ip defined as the angle in 
the xy-plane from the £-axis) and a polar (# defined as the 
angle form the z-axis) component so that the unit vectors 
denote a right handed system with e r A = e$. In Fig. U 
and [3] we plot the probability density functions for the radial 
respectively (^-component of the velocity field in linear scale; 
omitting the plot for the mean velocity as it does not show 
qualitatively different features than the mean y?- velocity. 

Often it is assumed that bulk velocities in dark mat- 
ter haloes are zero. Again, this is approximately true in the 
central part. The standard deviation of the local mean ra- 
dial velocity is about 4% of the shell averaged radial velocity 
dispersion. In the outer parts, not even the spherically av- 
eraged values are zero (see also Table [I]) and the local bulk 
velocities reach up to twice the values of the local velocity 
dispersion. At 400 kpc, the dispersion of the mean radial 



velocity exceeds 40% of the shell averaged radial velocity 
dispersion. 

At different radii one can see two populations of spheres: 
one that is infalling and another one that is moving out- 
wards. This is best seen along the x-axis but often the two 
populations are too much smeared out. Such a distribu- 
tion naturally arises from a cosmological infall pattern and 



has been seen before in numerical simulations (Diemand & 
|Kuhlen||2008| ). 

A similar picture emerges from the mean ip- and 
velocities with a comparable spread of the probability den- 
sity distributions at a given distance as for the mean ra- 
dial velocities. We find that the sub-population of spheres 
along the long axis have a narrower and more centrally 
peaked distribution in mean <p- and ^-velocities than the 
sub-populations along the short and intermediate axes indi- 
cating some degree of coherent tangential flow around the 
short and intermediate axes. 

The distribution of the subhalo-affected sample of 
spheres is always much broader than the total sample, and it 
contributes most to the extremes (both high and low) of the 
overall distribution. This indicates that the subhalo popula- 
tion as a whole is kinematically hotter than the background, 



in agreement with earlier studies by Diemand et al. (2004). 



Indeed we find subhaloes in the inner halo that have a speed 
of around 500 km s _1 w 3.5 cr r (8 kpc). 

In the inner region, we find that the subhalo-affected 
sample is skewed towards negative radial velocities (i.e. in- 
falling spheres) and positive (^-velocities. The skew in radial 
velocity is a signature of the tidal disruption process: we 
see these subhaloes on their last infall before they either are 
completely disrupted or they loose so much mass that they 
do not contribute much to the sphere averaged properties 
any more on their way out since they are now much more 
compact and less massive. 

All this conclusively shows that the traditional notion 
of a dark matter halo with no local bulk velocities is an inac- 
curate description for structures that form in cosmological 
simulations, especially in the outer parts of haloes. 



3.4 Local velocity dispersions 

In Fig. ^] [5] and [6] we plot the probability density functions of 
the local total, radial and (^-velocity dispersion. We omit the 
figure for the ^-velocity dispersion since it shows a similar 
behaviour as the ip- velocity dispersion. The local velocity 



dispersion is given by g\ — v\ — v'f, for k = r, ip or $ and 
&a + &b + ^"2 5 where cr a , (Jb respectively a c are the 
values of the dispersions in the eigencoordinate system of the 



-2 

Otot 



local dispersion ellipsoid (see section 3.6 for more details) 



We see that the local velocity dispersion has strong un- 
derlying variations. Whereas the central part only has varia- 
tions on the few per cent level with respect to the spherically 
averaged value (e.g. at 8 kpc we have a(a r )/a r ~ 0.05), one 
can find regions in the outer part of the halo that deviate by 
approximately an order of magnitude in velocity dispersion 
and at 400 kpc we have for example a(a r )/a r ~ 0.25. 

The total velocity dispersion shows the most compact 
probability density function when measured with respect to 
the spherically averaged value, i.e. o"(<rtot)/0tot is smaller 
than cr(afc)/o"fc for k = r, (p or $ and the distributions peak 
generally around the spherically averaged value. 
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Figure 2. Probability density functions of v r at different galactocentric distances r normalised by the spherically averaged radial velocity 
dispersion a r - We plot the same sub-samples as in Fig. ^ and use the same notation. 



Interestingly, we find that regions in the centre along the 
long axis are colder (i.e. a to t is smaller) than spheres along 
the intermediate or short axes. This trend is not observed 
in the outskirts of the halo. We also find that along the long 
axis, the radial velocity dispersions are higher and the tan- 
gential velocity dispersions are low in the inner region of the 
halo. This is due to the prolate shape and orientation of the 
local velocity dispersion ellipsoid and is discussed in more 
detail in section |3l3l In general, we find that the distributions 
of velocity dispersions for the two tangential velocity com- 
ponents are much broader along the short axis than along 
the long and intermediate axes. The subhalo- affected sample 



shows a less peaked distribution in velocity dispersions than 
the total sample. In the regions with enough subhaloes, the 
distribution is skewed towards lower dispersions for all ve- 
locity components. This is of course due to the lower internal 
velocity dispersion in subhaloes. 

It remains to be seen if these empirical trends of veloc- 
ity dispersion with orientation are universal for dark matter 
haloes or if they depend on the detailed hierarchical build-up 
history of every individual halo. Nevertheless, it is clear that 
the degree of chaotic motion can strongly vary locally and 
deviate by substantial factors from the spherically averaged 
value. 
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Figure 3. Probability density functions of at different galactocentric distances r normalised by the spherically averaged ^-velocity 
dispersion cr^. We plot the same sub-samples as in Fig. [l] and use the same notation. 



3.5 Local anisotropy parameter 

In Fig. [7| we plot the probability density function of the 
local anisotropy parameter /3 within the spheres. The local 
anisotropy parameter is given by 



15? 
2 al 



(8) 



where a\ — + o# is the tangential velocity dispersion 
squared and where we neglect the correlation term. 

This form of the anisotropy parameter is often used but 
has some problems since the following assumptions which 



are hidden in this expression are in general not fulfilled: i) 
v% — ii) v r = v^p — v$ — and hi) cov(^,^) = 0. 
These three conditions are only approximately fulfilled in 
the central part of a halo and certainly not in the outskirts. 
One should therefore see the expression for (3 as a definition 
for local anisotropy so that a comparison with previous work 
is possible. 

We find that the radial anisotropy of the velocity dis- 
persion tensor in our halo increases with radius (see also 
Table [T]), in agreement with previous studies (see e.g. 



Cole 



fe Lacey|1996||Colm et al.|2000| |Fukushige fc Makino| 2Q01 . 
Diemand et al.|2004| |Hansen & Moore|2006| ). In the central 
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Figure 4. Probability density functions of the total velocity dispersion <7tot at different galactocentric distances r normalised to the 
spherically averaged value a"tot- We plot the same sub-samples as in Fig. [I] and use the same notation. 



part we find that regions along the long axis are preferen- 
tially on radial orbits, whereas regions along the short axis 
are tendencially on tangential orbits. This effect disappears 
at larger radii. It can be understood as a direct consequence 
of the variation of the orientation of the local velocity dis- 
persion ellipsoid discussed in the next section [376] 

The probability densities peak in general close to the 
spherical average value except at 400 kpc where it is peaked 
towards highly radial anisotropics. This is especially true for 
the the sub-profiles of the y- and z-axis sample whereas the 
x-axis sample peaks around /3 = 0. With increasing distance 
also the spread in local anisotropy becomes larger and can 



vary from close to perfectly radial to close to perfectly tan- 
gential. By inspecting the subhalo- affected sample, we find 
that this sample is skewed towards tangential values of (3 at 
galactocentric distances with enough subhaloes. 



3.6 Local velocity dispersion ellipsoid 

In each sphere we also calculate the local velocity covariance 
tensor, also know as the velocity dispersion tensor, X) given 
by 

Sy EE (Vi - Vi){Vj -Vj) . (9) 
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Figure 5. Probability density functions of radial velocity dispersion ay at different galactocentric distances r normalised to the spherically 
averaged value ay. We plot the same sub-samples as in Fig. ^ and use the same notation. 



By diagonalising XI, the velocity dispersions in the eigenco- 
ordinate system are given by the square roots of the eigenval- 
ues of the velocity covariance tensor. These velocity disper- 
sions in the eigencoordinate system define the local velocity 
dispersion ellipsoid and are sorted by a a ^ &b ^ &c- We 
name the appropriate eigenvectors accordingly, e.g. e aa etc. 

First, we check the orientation of the local velocity dis- 
persion ellipsoid. For this purpose we calculate the angles 
between the eigenvectors of the velocity covariance tensor 
and the appropriate shape axis given by 

a a = arccos(e x • e aa ) (10) 



otb = arccos(e y • e ah ) (11) 
a c = arccos(e 2 • e<j c ) (12) 

In Fig. [8] respectively [5] we plot the probability density func- 
tion of a a and a&. We do not show the figure for a c since 
it shows a qualitatively similar behaviour as for ab- We also 
plot the curve, which is proportional to sin(afc), that corre- 
sponds to an isotropic distribution. 

We observe that in the inner regions the local veloc- 
ity dispersion ellipsoid is close to perfectly aligned with the 
shape ellipsoid since all angle probability density functions 
peak sharply at small angles. This effect is most prominent 
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Figure 6. Probability density functions of (^-velocity dispersion at different galactocentric distances r normalised to the spherically 
averaged value a^. We plot the same sub-samples as in Fig. fl] and use the same notation. 



for the major axis sub-sample, which has a very tight and 
sharp distribution around small angles ak for all k, and 
less so for the intermediate and short axes, which in some 
cases show a broader spread towards larger angles. The gen- 
eral alignment of the local velocity dispersion ellipsoid with 
the shape of the dark matter halo nicely explains the ob- 
served behaviour of the anisotropy parameter /5, for which 
we mainly found radial orbits along the x-axis, isotropic 
orbits along the intermediated y-axis and tangential orbits 
along the z-axis. The further out we go, the more isotropic 
the different angle distributions become. However, at nearly 
all distances from the Galaxy centre we find that smaller an- 



gles are slightly more probable than in a perfectly isotropic 
distribution. 

Interesting is the behaviour of the different sub-samples 
further out. The alignment with the shape seems to be best 
along the x-axis. Only in the outskirts at 400 kpc, this sub- 
sample distribution becomes more and more isotropic for all 
angles. The sub-samples along the other two axes show even 
in the inner part a broader distribution that does not fol- 
low the overall distribution. For example, the z-axis sample 
shows a distribution peaked around large angles for ab and 
a c for intermediate distances around 50 kpc and 100 kpc 
whereas the a a distribution is mildly peaked towards small 
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Figure 7. Probability density functions of the local anisotropy parameter (3 at different galactocentric distances r. We plot the same 
sub-samples as in Fig. ^ and use the same notation. 



angles at 50 kpc (not shown in figure). In other words: only 
the long axis seems to be slightly aligned with the long axis 
of the shape ellipsoid whereas the the intermediate and short 
axes are perpendicular to the corresponding shape axes. In 
general, also the subhalo sample shows a rather isotropic 
distribution although it also tends to be aligned with the 
shape ellipsoid in the inner part of the halo. 

Globally, we expect such an alignment from the tensor 



virial theorem (Binney & Tremaine 1987) and such a cor- 



relation has previously been found in cosmological iV-body 
simulations, e.g. in Allgood et al. (2006). If we calculate the 



velocity dispersion ellipsoid for all the particles in the shell, 



we also find a very tight alignment of this shell dispersion 
ellipsoid with the shape ellipsoid (see e.g. Fig. [^respectively 
|9| . Only in the outer parts at 400 kpc we get a deviation from 
that behaviour. But in principle, we do not expect that the 
tensor virial theorem holds locally. It seems that the align- 
ment of the local velocity dispersion ellipsoid with the shape 
ellipsoid only holds in relaxed and well mixed regions - the 
central region of the Via Lactea II halo. Hence, this local 
alignment might just be a numerical artefact since the cen- 
tral region is too relaxed due to numerical under-resolving if 
compared to the true degree of relaxation in reality. But an 
alignment of the local dispersion ellipsoid with the shape is 
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Figure 8. Probability density functions of the angle a a at different galactocentric distances r. We plot the same sub-samples as in Fig. 
[Hand use the same notation. 



also found in observations. For example, galaxies from the 
SAURON survey that were dynamically modelled with the 
Schwarzschild technique also show such a correlation ( |Cap- 
pellari et aL|2007| ). 

Furthermore, we calculate the shape of the local velocity 
dispersion ellipsoids which we measure with the triaxiality 
parameter < |Franx et al.|1991| l deBned by 



■si 



(13) 



Ellipsoids are called oblate ifO^T^l/3, t riax ial if 1/3 < 
T < 2/3 and prolate if 2/3 < T < 1. In Fig. [To] we plot the 



probability density function of the triaxiality parameter in 
linear scale. 

The total probability density function shifts from a peak 
in the oblate region (e.g. at 25 kpc) via a peak in the tri- 
axial region (e.g. at 100 kpc) to a more prolate shape at 
400 kpc. Only the innermost region at 8 kpc does not fol- 
low this overall trend. Interesting is again the distribution of 
the different sub-samples. The shape of the velocity disper- 
sion ellipsoids along the x-axis generally peak in the prolate 
region although the distribution becomes broader towards 
oblate shapes further out. The ?/-axis sample has in the in- 
ner region a rather oblate shape which drifts via triaxial 
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Figure 9. Probability density functions of the angle at different galactocentric distances r. We plot the same sub-samples as in Fig. 
[Hand use the same notation. 



(e.g. at 100 kpc) to a prolate distribution in the outskirts of 
the halo. Completely different is the behaviour of the z-axis 
sample: at 8 kpc it's distribution is peaked in the prolate re- 
gion then drifts via the triaxial region at 25 kpc to a oblate 
distribution at 50 kpc (not shown in the figure). Then it 
swings back via a rather triaxial shape to a prolate shape 
again at 400 kpc. The shape of the shell averaged velocity 
dispersion tensor does not fully follow the trend of the over- 
all distribution and is in general not close to the peak of the 
total probability density function. 

At the moment it is not clear what the underlying cause 
for these shape and orientation variations with galactocen- 



tric distance is and if these trends are universal. Further 
numerical investigations are necessary. 

3.7 Local velocity space 

We now turn to a closer look at the local velocity space 
structure. For this purpose, we present 3D visualizations of 
the velocity structure within spheres cut out along the y- 
axis, with radii given by r sp h(r) from Table [I] 

Fig. [XT] shows the positions and velocity vectors for ev- 
ery particle within these spheres (1st and 3rd row) for all 
six galactocentric distances we investigated. We also plot the 
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Figure 10. Probability density functions of the triaxiality parameter T at different galactocentric distances r. We plot the same sub- 
samples as in Fig. ^ and use the same notation. 



location of the particles in the local velocity space (2nd and 
4th row). The colour (and length in position space) encodes 
the magnitude of the velocity vector. The big white arrow 
in the centre of the position space plots points towards the 
galactic centre whereas the black cube in the centre of the 
local velocity space plots marks the origin. For the veloc- 
ity space plots, we split the velocity vector field in a radial- 
(a;-axis), ip- (2/- axis) and a ^-component (2-axis). 

These spheres have been selected to be free of subhaloes. 
Therefore, one would not expect to see clumpy structure in 
the phase space of these spheres. In the inner spheres no 
velocity space structure is apparent. The orientation of the 



velocity vectors appears random and the actual velocity dis- 
tribution can be well fit by a multivariate normal distribu- 
tion (the different velocity dispersion components are given 
in Table [I]). But further out, there is evidence for locally 
coherent motion, visible as groups of vectors with the same 
colour pointing in the same direction. This is even clearer 
in the velocity space plots, which exhibit an increasing de- 
gree of dumpiness the further out one goes. At 400 kpc, 
for example, there are only a bit more than a dozen clumps 
in velocity space and no smooth component at all. Obvi- 
ously a smooth multivariate normal distribution would not 
provide a good fit. This is evidence that the outer part of 
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Figure 11. Positions of particles within the selected spheres along the y-axis with velocity vector attached (1st and 3rd row) as well 
as the corresponding location of the particles in the local velocity space (2nd and 4th row). The colour in both cases is given by the 
magnitude of the velocity vector. 
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Figure 12. Central sphere of radius 50 kpc where we calculated the position space density (left) and the true phase space density (right) 
with EnBiD. Approximately 2000 peaks from the subhaloes within that radius are visible in both pictures but the contrast of subhaloes 
is higher in the phase space density picture. In the phase space picture are also large scale dark matter streams visible which were formed 
by tidal mass loss of infalling subhaloes. These streams are not visible in the traditional position space density picture. 



the halo is built up from a collection of large scale streams 
from the tidal disruption of infalling subhaloes. But so too 
is the smooth component in centre, in that it also likely con- 
sists of an overlap of many, many streams. The central limit 
theorem then guarantees that the resulting distribution in 
the centre closely resembles something like a multivariate 
normal distribution (|Helmi et aL]|2003| or a generalisation 
thereof ( |Hansen et al.|2006 ). 

In Fig. |12| we show the central region of Via Lactea 
II: the sphere has a radius of 50 kpc. We plot the position 
space density (left) and t he true phase space density ( right) 
calculated with EnBiCp] ( [Snarma fe Stemmetz||2006} . It is 
important here to calculate the true phase space density in 
6D and not the commonly used pseudo phase space density 
p/a 3 (where a is the velocity dispersion), since in the latter 
information in velocity space is lost by averaging over the 
particles in local position space, instead of calculating the 
density in local velocity space directly from the particles. 
We only show the top 5 to 6 orders of magnitude in position 
and phase space density, in order not to overload the two 
pictures. 

Many cold streams are clearly visible in the central re- 
gion. Although these cold streams only contribute a few per- 
cent to the local mass density, their velocity dispersion is just 
a few km s _1 , resulting in a very high phase space density for 
these particles. These streams are not visible in traditional 
density or density squared pictures and can only be revealed 
by visualising the true phase space density. Via Lactea II 
seems to be the first structure formation simulation with 
sufficient resolution to reveal these streams in phase space, 



http : //sourcef orge . net/pro j ects/enbid/ 



since previous lower resolution runs did not show such phase 
space features. Also subhaloes are better visible due to their 
higher contrast in phase space density. Approximately 2000 
peaks are seen in the central 50 kpc sphere phase space den- 
sity image. The high contrast makes the phase space density 
the ideal method of finding subhaloes. 

From our findings it is obvious that also the central 
region shows a huge amount of additional structure aside 
from the expected subhalo density peaks. A more detailed 
discussion about streams and their properties will follow in 
a future publication. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we present a study of local properties of the 
Via Lactea II dark matter halo. Commonly, characteristics 
of dark matter haloes are described by spherically symmetric 
profiles. Here, we show that locally, at a fixed galactocen- 
tric distance, these properties can vary by orders of magni- 
tude from their canonical, spherically averaged values. For 
example, while the density is smooth in the centre, in the 
outskirts the density range spans four orders of magnitude. 
This is due to the presence of both subhaloes and underdense 
regions or holes in the matter distribution. The widespread 
assumption that there are no bulk flows is only warranted 
in the central region, where we also find the local velocity 
dispersion ellipsoid to be aligned with the shape ellipsoid of 
the halo. Neither of these findings hold in the outer parts of 
the halo, where particles exhibit bulk motions and a more 
isotropic distribution of the velocity dispersion ellipsoid's 
orientation is found. The local velocity space structure can, 
in general, not be well described by a smooth multivariate 
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normal distribution, at least not in the outer parts of the 
halo. A qualitatively new feature in Via Lactea II is the de- 
tection of streams that are made visible only through their 
true 6D phase space density. 

Such variations of local properties is due to both the tri- 
axial shape of the dark matter halo and its generally clumpy 
structure in phase space. We find that the phase space struc- 
ture of a dark matter halo shows a significant departure from 
the canonical picture of a smooth background density pro- 
file with subhaloes in position space and multivariate normal 
distributions in velocity space. Spherically averaged quanti- 
ties do not adequately describe properties of the dark matter 
halo at a given galactocentric distance - especially not the 
wealth of structure we find locally - since a lot of information 
gets smeared out by spherically averaging. In general, dark 
matter haloes show a high degree of graininess - clumps in 
phase space - which will probably show up at even smaller 
scales that now seem to be smooth in future simulations 
with higher resolution. Therefore, a smooth and featureless 
distribution function does not accurately describe dark mat- 
ter haloes that form in a cosmological structure formation 
simulation. 

We find several correlations between the shape ellipsoid 
and the local velocity dispersion ellipsoid. It is unclear at 
the moment to what degree our findings are universal or 
how these correlations depend on the hierarchical build up 
history. 



Knebe & WieBner (2006) earlier estimated analytically 



the effect of a triaxial halo shape on the variance of the lo- 
cal density. Their values of the dispersion normalised to the 
spherically averaged value for a halo with a similar triax- 
iality correspond quite good with our values in the centre 
of the halo but we get much higher dispersion values (larger 
than the mean) in the outskirts of the halo. This difference is 
due to their assumption of a smooth triaxial density profile 
and we get a higher scatter due to the presence of subhaloes 
and underdense regions. 

For dark matter direct detection experiments it is essen- 
tial to know the local dark matter properties at 8 kpc. The 
spherically averaged value for the density in Via Lactea II is 
p(8 kpc) = 1.056 x 10" 2 M pc~ 3 = 0.4008 GeV c~ 2 cm" 3 , 
close to 0.3 GeV c~ 2 cm -3 which is often used as a canon- 
ical value in the literature ( Kamionkowski & Kinkhabwala 



1998 Particle Data Group 2008). But as we have shown 



here, this value can vary locally. One of the large uncer- 
tainty factors is the missing information about the orienta- 
tion of the disk with respect to the dark matter halo. There 
are different claims from observations and theory about a 
possible alignment of the angular momentum axis of the 



disk with the shape axes of the halo (see e.g. 


Navarro et al. 


|2004| |Bailin et al.|2005| |Sharma & Steinmetz 


2005, and ref- 



erences therein) so that a clear answer is not possible at 
the moment. But often it is claimed that the angular mo- 
mentum axis of the disk and the short axis (z) tend to be 
aligned (Baili^eT^Lj2005 ) which would mean that the disk 



would lie preferentially in the xy-p\&ne in our coordinate 
system. An additional problem is that we measure these lo- 
cal properties within a sphere of r sp h(8 kpc) = 500 pc radius 
but for dark matter detection experiments more a scale of 
1 AU = 4.848 x 10" 6 pc = 9.696 x 10" 9 r sph (8 kpc) is rel- 
evant and it is not clear what the local properties of the 
dark matter distribution on a 1 AU scale are or how one 



could reasonably extrapolate that over 8 orders of magni- 
tude - especially considering the highly non-linear numeri- 
cal effects mentioned above and in the appendix that affect 
the local phase space structure. Also the missing baryonic 
physics in Via Lactea II like adiabatic contraction, stellar 
disk and bulge, inspir ailing compact objects like black holes 
etc. can modify the central dark matter structure in either 
way. Therefore, it is still not clear what the detailed struc- 
ture of the dark matter locally is. 
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APPENDIX A: COMPARISONS 

In this section we investigate the influence of different defi- 
nitions of locality and numerical resolution on our local es- 
timates. 



Al Definition of locality 

Our definition of locality is somewhat arbitrary. If the re- 
gion is too big then local properties get smeared out by the 
averaging procedure and if the region is too small then sta- 
tistical fluctuations become important. Hence, we chose the 
size of our spheres so that they contain (9(10 3 ) particles as in 
similar previous work ( Moore et al.|200T Helmi et al.|2 002). 
Exactly: an average density sphere of radius 1 r sp h contains 
1356 particles. In order to see the effect of different choices 
for our local estimate, we repeated the exercise from section 
[3] with spheres of different radii of the following size: 1/4 r sp h, 
1/2 r sp h, and 2 r sp h, with a corresponding increase or de- 
crease in the number of particles sampled per sphere. The 
resulting Poisson scatter ranges from ~22% for the smallest 
spheres (1/4 r sp h) to ~1% for the largest spheres (2 r sp h) 
if we take the spherical averaged density as a reference. By 
only looking at the lowest density along the z-axis we get 
a range of ~41% for the smallest spheres to ~1.7% for the 
largest spheres. 

In Fig. |A1| we present the effect of the different defi- 
nitions of locality on the local density distribution and see 
that they show the expected result. The regions that contain 
8 times more particles are of course smoothed to a higher 
degree so that the resulting spread in the probability distri- 
bution function is a bit reduced. In a similar way the spread 
for the 8 times smaller spheres is increased. But the different 
definitions of locality mainly affect the rare outliers on the 
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Figure Al. Probability density functions of the local density p at different galactocentric distances r for different definitions of locality 
normalised by the standard spherically averaged value p given in Table [l] An average density sphere of radius 1 r sp h contains 1356 
particles. The 1/4 r sp h probability density function is significantly broadened by Poisson noise. 



tails of the distribution. By reducing to 64 times less par- 
ticles, the extremes of the distribution become much more 
populated, which is mostly due to the increased Poisson scat- 
ter and to a lesser extend due to increased graininess on the 
smaller probed scales. 

Any size of a sphere that contains at least a few hundred 
particles would be fine for a definition of locality, indicating 
that our choice is a rather conservative choice and the result- 
ing variances of the different distributions from the previous 
sections are lower limits. The probability density functions 
for the other characteristics show the same effects as func- 



tion of locality definition as the one discussed here for the 
density. 



A2 Influence of resolution 

We compare the results from Via Lactea II (VL2) to a 
medium resolution simulation (VL2m) of the same halo. In 
the medium resolution run, the particle mass was a factor 
64 higher than in the high resolution run. Therefore, in or- 
der to probe local properties with approximately the same 
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Figure A2. Probability density functions of the local density p for different resolutions of Via Lactea II normalised by the spherically 
averaged value p for the high resolution run given in Table [l] An average density sphere of radius 1 r sp h contains 1356 particles in VL2, 
and 21 particles in VL2m. This means the 1 r sp h measurements in the VL2m simulation are rather noisy. 



number of particles spheres with \/(34 — 4 times larger radii 
were used. 

In Fig. |A2| we show the probability density functions 
of the local density for the medium and high resolution run 
for different sizes of the spheres. In general the rarer peaks 
at the low and high end of the distribution are missing in 
the medium resolution run when one compares spheres with 
equal number of particles in both simulations (i.e. VL2 with 
1 r sp h and VL2m with 4 r sp h) since the low resolution run 
resolves less substructure and it is smoothed over a 4 times 
larger scale. When comparing spheres of equal physical size 



in both runs, then the low resolution run shows a broader 
distribution due to additional Poisson noise. 

The lack of high density peaks in the low resolution sim- 
ulation is due to numerical effects. Less subhaloes survive in 
the medium resolution run since the subhaloes are resolved 
with fewer particles and are more easily tidally disrupted. 
This leads to the effect that a larger region in the centre is 
smooth in lower resolution runs. The higher the resolution 
the smaller this apparently smooth region is. 

Additionally, the higher degree of numerical artefacts 
(e.g. heating, relaxation etc.) in the medium resolution run 
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smears out streams. For example, no dark matter streams 
are visible in the inner 50 kpc in phase space density maps of 
the medium resolution of Via Lactea II and the peak phase 
space densities of the subhaloes are much lower than the 
values from the high resolution run. 

It is not clear at the moment to what degree the dark 
matter streams in the high resolution run are broadened 
by artificial heating in the simulations. Of course there are 
also real dynamical heating sources, such as dark matter 
subhaloes or baryonic structures like a stellar disk, a bulge, 



or gas clouds (see e.g. Ibata et al. 2002 Siegal-Gaskins & 
|Valluri|2008) . 

Material that ends up in the central regions of dark 



matter haloes originates from early forming halos (Diemand 



et al. 2005) and was accreted into the main halo early on 
( Helmi et al.|[2002 ). One therefore expects a higher degree 



of phase mixing due to the early accretion and short time- 
scales in the centre. Nevertheless it is not clear if the smooth 
appearance of the central regions in Via Lactea II is entirely 
due to effcient phase mixing: The small, early progenitors 
that build up the central dark matter halo of Via Lactea II 
are under resolved and therefore too low numerical resolu- 
tion might appear as an efficient phase mixing. 



empty bins (or spheres in our case), i.e. k — 0. This results 
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where A = l/n is the expectation value of balls per bin. We 
can simplify this in the large n limit to 
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which we use to estimate the expected fraction of empty 
bins respectively spheres. 



APPENDIX B: BALLS- IN- BINS STATISTICS 

Let's try to solve the following problem. We have n bins 
and balls are thrown randomly into the bins. What is the 
probability p(k, I) of finding bins that contain k balls after 
throwing / balls? 

The probability of a ball hitting a bin is given by pu = 
l/n and the probability of missing a bin is given by p m = 
l—ph = I — l/n. Obviously, we have p(0, 0) = 1, p(k, 0) = 
for k > and p(k, I) = for k > I. After throwing / — 1 balls 
we can write 

p(k, I) = Pmp{k, I - 1) + p h p(k -1,1-1) . (Bl) 

This is a recursion formula that allows us now to construct 
the general form of p(k, I). By setting p{— 1, I) = for all I, 
we get for the non-zero values for the first few values of / 
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and so on. 


It becomes obvious that the 


general pattern is 


given by a 


binomial distribution 
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We are specifically interested in the case of getting 



